1 Overview

This script visualizes and analyzes the mean amplitude of response (meanbeta) to expected and unexpected events from three “domains” - psychology-action, psychology-environment, and physics. Our pre-registered confirmatory analyses center around the psychology-action and physics videos.

  • Part 1a: Domain and event effects 8 focal regions. This is the univariate confirmatory analysis for the project.
  • Part 1b: Checking robustness of results depending on ROI size (pre-registered as 100 voxels)
  • Part 2: Repeating the univariate analysis in a larger set of domain-general and domain-specific regions, with more conservative significance thresholds given the # of regions we’re exploring per category.
  • Part 3: Exploratory univariate analysis on VOE effect in action-environment scenarios.
  • Part 4: Similar to confirmatory analysis, except that we use the VOE events from runs 3-4 to select ROIs. Function of this analysis is to check for whether the confusing STS results we obtained is due to conceptual mismatch between the DOTS localizer and the psych-action events.
  • Part 5: We model betas for each scenario and event, rather than folding across scenarios / tasks. This allows us to test for effects of visual statistics (Part 5a), which were computed over each video individually. Also allows us to do per-task analyses (Part 5b).
  • Part 6: Compute and plot Dice coefficient between two focal regions (bilateral APC and L/RSMG).
  • Part 7: Region by region univariate effect size (analogue to confirmatory MVPA analysis that is rendered uninterpretable, given the lack of reliable VOE spatial patterns across runs)
knitr::opts_chunk$set(echo = TRUE)
options(scipen = 1, digits = 3)
library(pacman)
pacman::p_load(tidyverse,
               patchwork,
               papaja,
               ggforce,
               conflicted,
               here,
               cowplot,
               lme4,
               lmerTest,
               effects,
               effectsize,
               EMAtools,
               performance,
               parameters,
               BayesFactor,
               gt,
               insight,
               sjPlot,
               lsmeans,
               ggrepel,
               GGally,
               nptest,
               boot,
               boot.pval,
               car)

conflict_prefer("select", "dplyr")
## [conflicted] Will prefer dplyr::select over any other package
conflict_prefer("filter", "dplyr")
## [conflicted] Will prefer dplyr::filter over any other package
conflict_prefer("rename", "dplyr")
## [conflicted] Will prefer dplyr::rename over any other package
conflict_prefer("mutate", "dplyr")
## [conflicted] Will prefer dplyr::mutate over any other package
conflict_prefer("summarise", "dplyr")
## [conflicted] Will prefer dplyr::summarise over any other package
conflict_prefer("lmer", "lmerTest")
## [conflicted] Will prefer lmerTest::lmer over any other package
conflict_prefer("here", "here")
## [conflicted] Will prefer here::here over any other package
source(here("helper_funs.R"))

options(contrasts = c("contr.sum", "contr.poly")) 

theme_set(theme_cowplot(10))

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] car_3.0-8              boot.pval_0.4          boot_1.3-25           
##  [4] nptest_1.0-3           GGally_2.1.2           ggrepel_0.8.2         
##  [7] lsmeans_2.30-0         emmeans_1.4.8          sjPlot_2.8.10         
## [10] insight_0.15.0         gt_0.8.0               BayesFactor_0.9.12-4.4
## [13] coda_0.19-3            parameters_0.16.0      performance_0.8.0     
## [16] EMAtools_0.1.3         effectsize_0.6.0.1     effects_4.1-4         
## [19] carData_3.0-4          lmerTest_3.1-2         lme4_1.1-23           
## [22] Matrix_1.2-18          cowplot_1.0.0          here_1.0.1            
## [25] conflicted_1.1.0       ggforce_0.3.4          papaja_0.1.0.9997     
## [28] tinylabels_0.2.3       patchwork_1.1.2        forcats_0.5.0         
## [31] stringr_1.5.0          dplyr_1.0.9            purrr_0.3.4           
## [34] readr_1.3.1            tidyr_1.2.0            tibble_3.1.7          
## [37] ggplot2_3.4.1          tidyverse_1.3.0        pacman_0.5.1          
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1        backports_1.4.1     plyr_1.8.6         
##   [4] splines_4.0.2       TH.data_1.0-10      digest_0.6.31      
##   [7] htmltools_0.5.4     fansi_1.0.3         magrittr_2.0.3     
##  [10] memoise_2.0.1       openxlsx_4.1.5      modelr_0.1.8       
##  [13] sandwich_2.5-1      colorspace_2.0-3    blob_1.2.1         
##  [16] rvest_0.3.6         mitools_2.4         rbibutils_2.2.8    
##  [19] haven_2.4.3         xfun_0.36           crayon_1.5.1       
##  [22] jsonlite_1.8.4      survival_3.1-12     zoo_1.8-8          
##  [25] glue_1.6.2          polyclip_1.10-0     gtable_0.3.0       
##  [28] MatrixModels_0.4-1  sjstats_0.18.0      sjmisc_2.8.5       
##  [31] abind_1.4-5         scales_1.2.0        mvtnorm_1.1-1      
##  [34] DBI_1.1.2           ggeffects_0.16.0    Rcpp_1.0.8         
##  [37] xtable_1.8-4        foreign_0.8-80      survey_4.0         
##  [40] datawizard_0.2.3    httr_1.4.2          RColorBrewer_1.1-3 
##  [43] ellipsis_0.3.2      pkgconfig_2.0.3     reshape_0.8.8      
##  [46] farver_2.1.0        nnet_7.3-14         sass_0.4.4         
##  [49] dbplyr_1.4.4        utf8_1.2.2          tidyselect_1.1.2   
##  [52] rlang_1.0.6         munsell_0.5.0       cellranger_1.1.0   
##  [55] tools_4.0.2         cachem_1.0.6        cli_3.5.0          
##  [58] generics_0.1.2      sjlabelled_1.1.7    broom_0.8.0        
##  [61] evaluate_0.19       fastmap_1.1.0       yaml_2.3.6         
##  [64] knitr_1.41          fs_1.5.2            zip_2.2.0          
##  [67] pbapply_1.6-0       nlme_3.1-148        xml2_1.3.2         
##  [70] compiler_4.0.2      rstudioapi_0.13     curl_4.3           
##  [73] reprex_0.3.0        statmod_1.4.34      tweenr_1.0.2       
##  [76] bslib_0.4.2         stringi_1.7.8       lattice_0.20-41    
##  [79] nloptr_1.2.2.2      vctrs_0.5.1         pillar_1.7.0       
##  [82] lifecycle_1.0.3     Rdpack_2.3.1        jquerylib_0.1.4    
##  [85] estimability_1.3    data.table_1.13.0   R6_2.5.1           
##  [88] rio_0.5.16          codetools_0.2-16    MASS_7.3-51.6      
##  [91] assertthat_0.2.1    rprojroot_2.0.3     withr_2.5.0        
##  [94] multcomp_1.4-13     bayestestR_0.11.5   hms_0.5.3          
##  [97] DataCombine_0.2.21  grid_4.0.2          minqa_1.2.4        
## [100] rmarkdown_2.19.2    numDeriv_2016.8-1.1 lubridate_1.7.9
regioninfo <- read.csv(here("input_data/manyregions_info.csv")) 
DT::datatable(regioninfo %>% dplyr::arrange(ROI_category, focal_region), options = list(scrollX = TRUE))
domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")
univariate_data <-
  read_rds(here("outputs/univariate_data/study2_univariate_data_allregions_alltopN.Rds"))

univariate_data$event <- relevel(univariate_data$event, ref = "unexp") # for some reason this is required for the model coefficients to reflect the intuitive def that unexpected > expected is a positive coefficient

univariate_summary_domain <-
  readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain.Rds"))

univariate_summary_domain_runs12 <-
  readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_runs12.Rds"))

univariate_summary_domain_allruns <-
  readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_allruns.Rds"))

2 PART 1a: Focal regions - confirmatory analysis

focal_data <- univariate_data %>%
  filter(focal_region == 1)

focal_summary_domain <- univariate_summary_domain %>%
  filter(focal_region == 1)

focal_summary_domain_runs12 <- univariate_summary_domain_runs12 %>%
    filter(focal_region == 1)

# focal_summary_task_runs12 <- univariate_summary_task_runs12 %>%
    # filter(focal_region == 1)

regions <- levels(as.factor(focal_data$ROI_name_final))

2.1 Plots

First, here are plots of betas, run by run, to physics and psychology events (familiarization, expected, unexpected).

plot_univar_runbyrun("superiortemporal_L", -1, 4) | plot_univar_runbyrun("superiortemporal_R", -1, 4)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("supramarginal_L", -1, 4) | plot_univar_runbyrun("supramarginal_R", -1, 4)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("antParietal_bilateral", -1, 6) | plot_univar_runbyrun("precentral_inferiorfrontal_R", -1, 6)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("V1_bilateral", -1, 11) | plot_univar_runbyrun("MT_bilateral", -1, 11)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

2.1.1 By domain (input from domain / event-level betas)

These plots go into Figure 3 of the main text.

# DS_runs12_plots <- plot_univar_runs12("superiortemporal", 0, 4) |
# plot_univar_runs12("supramarginal", 0, 4) 
# 
# DG_runs12_plots <- plot_univar_runs12("MD", 0, 6) | plot_univar_runs12("early_visual", 0, 11)
# 
# DS_runs12_plots / DG_runs12_plots

(plot_univar_runs12("superiortemporal_L", 0, 4) | plot_univar_runs12("superiortemporal_R", 0, 4) )
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

(plot_univar_runs12("supramarginal_L", 0, 4) | plot_univar_runs12("supramarginal_R", 0, 4))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

(plot_univar_runs12("antParietal_bilateral", 0, 6) | plot_univar_runs12("precentral_inferiorfrontal_R", 0, 6))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_runs12("V1_bilateral", 0, 11) | plot_univar_runs12("MT_bilateral", 0, 11)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

2.1.2 Action-env (input from domain / event-level betas)

summary_domain_allruns <- read_rds(here("outputs/descriptive_summaries/study2_univariate_summary_domain_allruns.Rds"))

Univariate response to psych_env events, over all runs, for all focal regions.

STS <- plot_both_allruns("superiortemporal_", -1, 10) 
SMG <- plot_both_allruns("supramarginal_", -1, 10)

V1 <- plot_both_allruns("V1_bilateral", -1, 10) 
MT <- plot_both_allruns("MT_bilateral", -1, 10)

APC <- plot_both_allruns("antParietal_bilateral", -1, 10)
RFC <- plot_both_allruns("precentral_inferiorfrontal_R", -1, 10)

STS | SMG | V1 | MT | APC | RFC + plot_layout(widths = c(2, 2, 1, 1, 1, 1))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

2.1.3 By task (input from single betas)

study2_univariate_summary_task_singlebetas <-  readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_task_singlebetas.Rds")) %>%
  filter(focal_region == "1")
regions <- levels(as.factor(study2_univariate_summary_task_singlebetas$ROI_name_final))

plot_univar <- function(region, ymin, ymax) {
  plotobject <- ggplot(data = study2_univariate_summary_task_singlebetas %>% 
                         filter(str_detect(ROI_name_final, region)),
           aes(x = event, y = meanbeta, fill = task)) +
      geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
      geom_errorbar(
        aes(ymin = meanbeta - se, ymax = meanbeta + se),
        position = position_dodge(width = .9),
        width = .2,
        colour = "black"
      ) +
      theme_cowplot(10) +
      facet_wrap( ~ ROI_name_final + factor(task, c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint")) + extracted_run_number, nrow = 2) +
      # scale_fill_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
      ylab("Average beta (amplitude)") +
      xlab("Event") +
      theme(axis.text.x = element_text(
        angle = 90,
        vjust = 0.5,
        hjust = 1
    ),
    legend.position = "none") +
    coord_cartesian(ylim = c(ymin, ymax)) +
    ggtitle(paste0("ROI:", region))
  
  plotobject
}
plot_univar("superiortemporal_L", -2, 3)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("superiortemporal_R", -2, 4)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("supramarginal_L", -2, 4)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("supramarginal_R", -2, 4)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("V1_bilateral", -1, 11)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("MT_bilateral", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("antParietal_bilateral", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.

plot_univar("precentral_inferiorfrontal_R", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.

2.2 Check whether we should restrict to first 2 runs

Following the preregistration plan, we check whether we are able to take the second level cope values (averaging across 4 runs) as the main dependent measure, or whether we need to restrict the analysis to just the first 2 runs.

First we fit a model including an interaction between run number and event.

focal_data_checking <- focal_data %>%
  filter(event != "fam",
         top_n_voxels == "100",
         domain != "both",
         str_length(extracted_run_number) == 4) %>%
  filter(ROI_name_final != "V1_bilateral",
         ROI_name_final != "MT_bilateral")

focal_data_checking$fixation_position <- as.factor(focal_data_checking$fixation_position)

focal_data_checking$event <- relevel(focal_data_checking$event, ref = "unexp")


checkmodel <- lmer(
  data = focal_data_checking,
  formula = meanbeta ~ extracted_run_number * event + (1|subjectID) + (1|ROI_name_final)
)

checkmodel %>%
 apa_print() %>%
 apa_table()
(#tab:unnamed-chunk-14)
**
Term \(\hat{\beta}\) 95% CI \(t\) \(\mathit{df}\) \(p\)
Intercept 2.20 [1.25, 3.16] 4.51 6.96 .003
Extracted run number1 0.41 [0.27, 0.55] 5.83 3,004.27 < .001
Extracted run number2 0.15 [0.02, 0.29] 2.18 3,004.27 .030
Extracted run number3 -0.15 [-0.29, -0.01] -2.06 3,006.55 .039
Event1 0.09 [0.01, 0.17] 2.27 3,003.97 .023
Extracted run number1 \(\times\) Event1 0.06 [-0.08, 0.19] 0.78 3,003.97 .433
Extracted run number2 \(\times\) Event1 0.07 [-0.07, 0.21] 1.00 3,003.97 .318
Extracted run number3 \(\times\) Event1 -0.07 [-0.21, 0.07] -0.92 3,003.97 .356
summary(checkmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ extracted_run_number * event + (1 | subjectID) + (1 |  
##     ROI_name_final)
##    Data: focal_data_checking
## 
## REML criterion at convergence: 13783
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -7.973 -0.603 -0.053  0.550  5.440 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  subjectID      (Intercept) 1.19     1.09    
##  ROI_name_final (Intercept) 1.20     1.10    
##  Residual                   5.14     2.27    
## Number of obs: 3048, groups:  subjectID, 32; ROI_name_final, 6
## 
## Fixed effects:
##                               Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                     2.2045     0.4890    6.9622    4.51   0.0028
## extracted_run_number1           0.4138     0.0710 3004.2738    5.83  6.1e-09
## extracted_run_number2           0.1544     0.0710 3004.2738    2.18   0.0297
## extracted_run_number3          -0.1484     0.0720 3006.5502   -2.06   0.0393
## event1                          0.0934     0.0411 3003.9740    2.27   0.0230
## extracted_run_number1:event1    0.0556     0.0709 3003.9740    0.78   0.4332
## extracted_run_number2:event1    0.0709     0.0709 3003.9740    1.00   0.3180
## extracted_run_number3:event1   -0.0662     0.0717 3003.9740   -0.92   0.3557
##                                 
## (Intercept)                  ** 
## extracted_run_number1        ***
## extracted_run_number2        *  
## extracted_run_number3        *  
## event1                       *  
## extracted_run_number1:event1    
## extracted_run_number2:event1    
## extracted_run_number3:event1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ext__1 ext__2 ext__3 event1 e__1:1 e__2:1
## extrctd_r_1 -0.001                                          
## extrctd_r_2 -0.001 -0.329                                   
## extrctd_r_3  0.002 -0.338 -0.338                            
## event1       0.000  0.000  0.000  0.000                     
## extrct__1:1  0.000  0.000  0.000  0.000 -0.005              
## extrct__2:1  0.000  0.000  0.000  0.000 -0.005 -0.330       
## extrct__3:1  0.000  0.000  0.000  0.000  0.014 -0.337 -0.337
lsmeans(checkmodel, pairwise ~ event | extracted_run_number)$contrasts
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3048' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3048)' or larger];
## but be warned that this may result in large computation time and memory use.
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3048' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3048)' or larger];
## but be warned that this may result in large computation time and memory use.
## extracted_run_number = run1:
##  contrast    estimate    SE  df z.ratio p.value
##  unexp - exp    0.298 0.164 Inf 1.821   0.0690 
## 
## extracted_run_number = run2:
##  contrast    estimate    SE  df z.ratio p.value
##  unexp - exp    0.328 0.164 Inf 2.008   0.0450 
## 
## extracted_run_number = run3:
##  contrast    estimate    SE  df z.ratio p.value
##  unexp - exp    0.054 0.166 Inf 0.327   0.7440 
## 
## extracted_run_number = run4:
##  contrast    estimate    SE  df z.ratio p.value
##  unexp - exp    0.066 0.164 Inf 0.405   0.6850 
## 
## Degrees-of-freedom method: asymptotic
model_habituation_results <- cbind(gen.m(checkmodel), gen.ci(checkmodel)[3:10,])
## Computing profile confidence intervals ...
model_habituation_results_byrun <- lsmeans(checkmodel, pairwise ~ event | extracted_run_number, pbkrtest.limit = 3048)$contrasts %>% as.data.frame()

We then compare this model to a simpler model without this interaction.

checkmodel_simpler <- lmer(
  data = focal_data_checking,
  formula = meanbeta ~ event + extracted_run_number + (1|subjectID) + (1|ROI_name_final)
)

summary(checkmodel_simpler)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event + extracted_run_number + (1 | subjectID) + (1 |  
##     ROI_name_final)
##    Data: focal_data_checking
## 
## REML criterion at convergence: 13774
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -7.947 -0.602 -0.057  0.553  5.414 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  subjectID      (Intercept) 1.19     1.09    
##  ROI_name_final (Intercept) 1.20     1.10    
##  Residual                   5.14     2.27    
## Number of obs: 3048, groups:  subjectID, 32; ROI_name_final, 6
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)              2.2045     0.4890    6.9620    4.51   0.0028 ** 
## event1                   0.0939     0.0411 3006.9741    2.29   0.0222 *  
## extracted_run_number1    0.4138     0.0710 3007.2741    5.83  6.1e-09 ***
## extracted_run_number2    0.1544     0.0710 3007.2741    2.18   0.0297 *  
## extracted_run_number3   -0.1484     0.0720 3009.5523   -2.06   0.0393 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 ext__1 ext__2
## event1       0.000                     
## extrctd_r_1 -0.001  0.000              
## extrctd_r_2 -0.001  0.000 -0.329       
## extrctd_r_3  0.002  0.000 -0.338 -0.338
# checkmodel_simpler %>%
#   apa_print() %>%
#   apa_table()

anova(checkmodel, checkmodel_simpler)
## refitting model(s) with ML (instead of REML)
knitr::kable(model_habituation_results)
B SE df t p CI_95_lower CI_95_upper
(Intercept) 2.204 0.489 6.96 4.508 0.003 2.208 2.323
extracted_run_number1 0.414 0.071 3004.27 5.831 0.000 1.183 3.226
extracted_run_number2 0.154 0.071 3004.27 2.175 0.030 0.275 0.553
extracted_run_number3 -0.148 0.072 3006.55 -2.062 0.039 0.015 0.293
event1 0.093 0.041 3003.97 2.274 0.023 -0.289 -0.007
extracted_run_number1:event1 0.056 0.071 3003.97 0.784 0.433 0.013 0.174
extracted_run_number2:event1 0.071 0.071 3003.97 0.999 0.318 -0.083 0.195
extracted_run_number3:event1 -0.066 0.072 3003.97 -0.924 0.356 -0.068 0.210
knitr::kable(model_habituation_results_byrun)
contrast extracted_run_number estimate SE df t.ratio p.value
unexp - exp run1 0.298 0.164 3004 1.821 0.069
unexp - exp run2 0.328 0.164 3004 2.008 0.045
unexp - exp run3 0.054 0.166 3004 0.327 0.744
unexp - exp run4 0.066 0.164 3004 0.405 0.685
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_perrun.Rds"))

This is a grey area. On the one hand we do not observe a significant run x event interaction. On the other hand, the primary VOE effect (unexp > expected) is indeed variable across runs, and stronger in runs 1 and 2 than runs 3 and 4, when considering psychology and physics events). So on balance, we have the grounds for focusing the analysis on runs 1-2.

Here we fit a model including main effects of domain and event (domain-general and domain-specific regions) and also their interaction (domain-specific regions only), and organize and print the outputs.

2.3 Model Tables

data_for_modeling <- univariate_data %>%
  filter(event != "fam",
         # top_n_voxels == "100",
         !str_detect(extracted_run_number, "fold"),
         !str_detect(extracted_run_number, "all_runs"),
         domain != "both")

ROIs <- levels(as.factor(data_for_modeling$ROI_name_final))
runs <- c("run1", "run2", "run3", "run4", "all_runs", "runs12")
top_n_voxels <- unique(data_for_modeling$top_n_voxels)
# top_n_voxels <- 100

effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)

rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "domain1", "event:domain")
modelsummaries_init <- data.frame()


for (ROI in ROIs) {
  curROI <- ROI
  
  for (topN in top_n_voxels) {
    cur_top_voxels_N <- topN
    
    for (run in runs) {
      curRun = run
      
      if (str_length(curRun) == 4) { # run1, run2, etc
        ROIdata <- data_for_modeling %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN,
                 extracted_run_number == curRun)
        
      } else if (curRun == "all_runs") {
        ROIdata <- data_for_modeling %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN,
                 str_length(extracted_run_number) == 4)
        
      } else if (curRun == "runs12") {
        ROIdata <- data_for_modeling %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN) %>%
          filter(extracted_run_number == "run1" | extracted_run_number == "run2")
      }
      
      # print("Current ROI: ", curROI)
      if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
        ROI_category <- "specific"
      } else {
        ROI_category <- "general"
      }
      
      # had to remove run random int (convergence))
      model <- lmer(data = ROIdata,
                    formula = meanbeta ~ event * domain + (1|subjectID))
      model_just_maineffects <-
        update(model, formula = ~ . - event:domain)  # Without interaction term
      model_just_domain <-
        update(model_just_maineffects, formula = ~ . - event)
      model_just_event <-
        update(model_just_maineffects, formula = ~ . - domain)
      
      BF_BIC_interaction <-
        data.frame(BF = exp((
          BIC(model_just_maineffects) - BIC(model)
        ) / 2),
        BF_interpretation = interpret_bf(exp((
          BIC(model_just_maineffects) - BIC(model)
        ) / 2)))  # BICs to Bayes factor
      
      BF_BIC_event <-
        data.frame(BF = exp((
          BIC(model_just_domain) - BIC(model_just_maineffects)
        ) / 2),
        BF_interpretation = interpret_bf(exp((
          BIC(model_just_domain) - BIC(model_just_maineffects)
        ) / 2)))
      
      BF_BIC_domain <-
        data.frame(BF = exp((
          BIC(model_just_event) - BIC(model_just_maineffects)
        ) / 2),
        BF_interpretation = interpret_bf(exp((
          BIC(model_just_event) - BIC(model_just_maineffects)
        ) / 2)))
      
      BF_initial <- rbind(BF_BIC_event,
                          BF_BIC_domain,
                          BF_BIC_interaction)
      
      rownames(BF_initial) <- rownames_no_intercept
      BFs <- rbind(BF_intercept_row,
                   BF_initial)
      
      effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
        select(d)
      effect_sizes <-
        rbind(effect_sizes_intercept_row, effect_sizes_initial)
      
      ROI_focal <- ROIdata$focal_region[1]

      modelsummary <-
        cbind(
          cur_top_voxels_N,
          curRun,
          ROI_category,
          ROI_focal,
          curROI,
          rownames(summary(model)$coefficients),
          summary(model)$coefficients,
          confint.merMod(model)[3:6, ],
          effect_sizes,
          BFs
        ) %>%
        as.data.frame()
      
      colnames(modelsummary) <-
        c(
          "top_n_voxels",
          "extracted_run_number",
          "domain",
          "focal_region",
          "ROI",
          "effect",
          "B",
          "SE",
          "df",
          "t",
          "p",
          "CI_95_lower",
          "CI_95_upper",
          "d",
          "BF",
          "BF_interpretation"
        )
      
      modelsummaries_init <-
        rbind(modelsummaries_init, modelsummary)
      cat(
        tab_model(
          model,
          show.std = TRUE,
          show.se = TRUE,
          show.stat = TRUE,
          auto.label = TRUE,
          title = paste0("Confirmatory Univariate ROI results in: ", curROI, " ", curRun, " top voxels = ", cur_top_voxels_N)
        )$knitr,
        "  \n  \n"
      )
    }
  }
}
modelsummaries <- modelsummaries_init %>%
  mutate_at(c(7:15), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

modelsummaries$effect <-
  as.factor(modelsummaries$effect)

levels(modelsummaries$effect) <-
  c("intercept", "domain", "event", "event:domain")
rownames(modelsummaries) <- NULL

# write_rds(modelsummaries, here("outputs/univariate_results/univar_results_alltopN_allregions.Rds"))
modelsummaries <- read_rds(here("outputs/univariate_results/univar_results_alltopN_allregions.Rds"))
DT::datatable(modelsummaries %>% filter(focal_region == 1, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect), options = list(scrollX = TRUE))

Below we check, for each domain-specific region (left and right STS and SMG), whether a model including an interaction between domain and event fits better than a model without the interaction.

focal_data_100_runs12 <- focal_data %>% 
        filter(top_n_voxels == "100",
               domain != "both", 
               event != "fam") %>%
        filter(extracted_run_number == "run1" | extracted_run_number == "run2")
LSMG_model <- lmer(
      data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
      formula = meanbeta ~ event * domain + (1|subjectID)
    )

plot(allEffects(LSMG_model))

# check_model(LSMG_model)
check_outliers(LSMG_model)
## OK: No outliers detected.
summary(LSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | subjectID)
##    Data: focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L")
## 
## REML criterion at convergence: 948
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4287 -0.5510  0.0345  0.5874  2.6185 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 4.02     2.00    
##  Residual              1.58     1.26    
## Number of obs: 256, groups:  subjectID, 32
## 
## Fixed effects:
##                Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)      2.3899     0.3630  31.0000    6.58  2.4e-07 ***
## event1           0.0672     0.0787 221.0000    0.85   0.3942    
## domain1          0.2580     0.0787 221.0000    3.28   0.0012 ** 
## event1:domain1   0.1877     0.0787 221.0000    2.39   0.0178 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1
## event1      0.000               
## domain1     0.000  0.000        
## event1:dmn1 0.000  0.000  0.000
LSMG_results_bydomain_top100 <-
  lsmeans(LSMG_model, pairwise ~ event |
            domain)$contrasts %>% as.data.frame()

knitr::kable(LSMG_results_bydomain_top100)
contrast domain estimate SE df t.ratio p.value
unexp - exp physics 0.510 0.222 221 2.29 0.023
unexp - exp psychology -0.241 0.222 221 -1.08 0.280
# write_rds(LSMG_results_bydomain_top100, path = here("outputs/univariate_results/LSMG_results_bydomain_top100.Rds"))


LSMG_model_simpler <- lmer(
      data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
      formula = meanbeta ~ event + domain + (1|subjectID)
    )
# check_model(LSMG_model_simpler)
check_outliers(LSMG_model_simpler)
## OK: No outliers detected.
LSMG_model_comparison <- anova(LSMG_model_simpler, LSMG_model) %>%
  as.data.frame()
## refitting model(s) with ML (instead of REML)
# write_rds(LSMG_model_comparison, path = here("outputs/univariate_results/LSMG_model_comparison_top100.Rds"))
RSMG_model <- lmer(
      data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_R"),
      formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
    )
summary(RSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |  
##     subjectID)
##    Data: focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_R")
## 
## REML criterion at convergence: 996
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.437 -0.585  0.036  0.515  3.808 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 6.484    2.546   
##  extracted_run_number (Intercept) 0.098    0.313   
##  Residual                         1.829    1.353   
## Number of obs: 256, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)      2.5420     0.5087  14.6472    5.00  0.00017 ***
## event1           0.1613     0.0845 219.9999    1.91  0.05767 .  
## domain1          0.5037     0.0845 219.9999    5.96    1e-08 ***
## event1:domain1   0.0566     0.0845 219.9999    0.67  0.50410    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1
## event1      0.000               
## domain1     0.000  0.000        
## event1:dmn1 0.000  0.000  0.000
# 
# RSMG_model %>%
#   apa_print() %>%
#   apa_table()
# check_model(RSMG_model)
check_outliers(RSMG_model)
## OK: No outliers detected.
lsmeans(RSMG_model, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.436 0.239 220 1.822   0.0700 
## 
## domain = psychology:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.209 0.239 220 0.876   0.3820 
## 
## Degrees-of-freedom method: kenward-roger
RSTS_model <- lmer(
      data = focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_R"),
      formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
    )
summary(RSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |  
##     subjectID)
##    Data: 
## focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_R")
## 
## REML criterion at convergence: 960
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4296 -0.5228  0.0122  0.5633  3.0469 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 2.6191   1.618   
##  extracted_run_number (Intercept) 0.0135   0.116   
##  Residual                         1.7596   1.326   
## Number of obs: 256, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)      2.0199     0.3090  20.7114    6.54  1.9e-06 ***
## event1           0.1234     0.0829 220.0000    1.49   0.1381    
## domain1         -0.2421     0.0829 220.0000   -2.92   0.0039 ** 
## event1:domain1   0.0539     0.0829 220.0000    0.65   0.5160    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1
## event1      0.000               
## domain1     0.000  0.000        
## event1:dmn1 0.000  0.000  0.000
# RSTS_model %>%
#   apa_print() %>%
#   apa_table()

# check_model(RSTS_model)
check_outliers(RSTS_model)
## OK: No outliers detected.
lsmeans(RSTS_model, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.355 0.234 220 1.512   0.1320 
## 
## domain = psychology:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.139 0.234 220 0.592   0.5540 
## 
## Degrees-of-freedom method: kenward-roger
LSTS_model <- lmer(
      data = focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_L"),
      formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
    )
summary(LSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |  
##     subjectID)
##    Data: 
## focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_L")
## 
## REML criterion at convergence: 1041
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.539 -0.593  0.035  0.503  3.454 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 4.337    2.083   
##  extracted_run_number (Intercept) 0.142    0.376   
##  Residual                         2.362    1.537   
## Number of obs: 256, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)      0.6181     0.4644   6.5575    1.33    0.228   
## event1           0.0624     0.0960 220.0000    0.65    0.517   
## domain1         -0.2998     0.0960 220.0000   -3.12    0.002 **
## event1:domain1  -0.0658     0.0960 220.0000   -0.69    0.494   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1
## event1      0.000               
## domain1     0.000  0.000        
## event1:dmn1 0.000  0.000  0.000
# LSTS_model %>%
#   apa_print() %>%
#   apa_table()

# check_model(LSTS_model)
check_outliers(LSTS_model)
## Warning: 1 outliers detected (cases 52).
outliers_list <- check_outliers(LSTS_model)
insight::get_data(LSTS_model)[outliers_list, ]
LSTS_model_nooutliers <- lmer(
      data = insight::get_data(LSTS_model) %>%
        filter(meanbeta != -11.6),
      formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
    )
# summary(LSTS_model_nooutliers)
LSTS_model_nooutliers %>%
  apa_print() %>%
  apa_table()
(#tab:unnamed-chunk-27)
**
Term \(\hat{\beta}\) 95% CI \(t\) \(\mathit{df}\) \(p\)
Intercept 0.62 [-0.29, 1.53] 1.33 6.56 .228
Event1 0.06 [-0.13, 0.25] 0.65 220.00 .517
Domain1 -0.30 [-0.49, -0.11] -3.12 220.00 .002
Event1 \(\times\) Domain1 -0.07 [-0.25, 0.12] -0.69 220.00 .494
lsmeans(LSTS_model_nooutliers, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp  -0.0069 0.272 220 -0.026  0.9800 
## 
## domain = psychology:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp   0.2564 0.272 220  0.944  0.3460 
## 
## Degrees-of-freedom method: kenward-roger

3 PART 1b: Sensitivity to voxel N

focal_modelsummaries <- modelsummaries %>%
  filter(focal_region == 1)

focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "runs12", "all_runs"))

ggplot(focal_modelsummaries %>%
         filter(str_length(extracted_run_number) == 4,
                effect != "intercept"), aes(extracted_run_number, d, colour = top_n_voxels)) +
  geom_point() +
  geom_line(aes(group = top_n_voxels)) +
  ylab("Cohen's D") +
  facet_wrap(~ effect + factor(ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_bilateral", "MT_bilateral")), nrow = 3) +
  # coord_cartesian(ylim = c(-.75, .75)) +
  geom_hline(yintercept = 0) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

4 PART 2: All regions (exploratory)

Then we repeat the same modeling procedure on these regions. But this time, we fit an interaction between event and domain in all regions, because some MD regions appear to have a preference for physical events, and it is an open question, if this is true, whether they too could encode physical prediction error.

DT::datatable(modelsummaries %>% filter(focal_region == 1, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect))

This set of regions includes left and right V1, MT, insula, and IFG (rather than bilateral). The bilateral versions of these ROIs are not included in this portion of the analysis.

DT::datatable(modelsummaries %>% filter(focal_region == 0, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect), options = list(scrollX = TRUE))

5 PART 3: Action-env scenarios

5.1 Check for run X event effects

focal_data_both_100 <- focal_data %>%
  filter(event != "fam",
         top_n_voxels == "100",
         domain == "both")
options(contrasts = c("contr.sum", "contr.poly"))

checkdata_both <- focal_data_both_100 %>%
  filter(ROI_name_final != "V1_bilateral",
         ROI_name_final != "MT_bilateral",
         str_length(extracted_run_number) == 4)

checkdata_both$fixation_position <- as.factor(checkdata_both$fixation_position)

checkdata_both$event <- relevel(checkdata_both$event, ref = "unexp")

checkmodel_bothdomains <- lmer(
  data = checkdata_both,
  formula = meanbeta ~ event * extracted_run_number + (1|subjectID) + (1|ROI_name_final)
)

summary(checkmodel_bothdomains)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * extracted_run_number + (1 | subjectID) + (1 |  
##     ROI_name_final)
##    Data: checkdata_both
## 
## REML criterion at convergence: 7354
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -8.428 -0.547 -0.028  0.538  5.262 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  subjectID      (Intercept) 1.11     1.05    
##  ROI_name_final (Intercept) 1.34     1.16    
##  Residual                   6.83     2.61    
## Number of obs: 1524, groups:  subjectID, 32; ROI_name_final, 6
## 
## Fixed effects:
##                               Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                     2.0973     0.5121    6.6088    4.10   0.0052
## event1                          0.2175     0.0669 1479.8794    3.25   0.0012
## extracted_run_number1           0.5830     0.1157 1480.2519    5.04  5.2e-07
## extracted_run_number2           0.1838     0.1157 1480.2519    1.59   0.1123
## extracted_run_number3          -0.4030     0.1172 1483.0707   -3.44   0.0006
## event1:extracted_run_number1   -0.1024     0.1156 1479.8794   -0.89   0.3760
## event1:extracted_run_number2   -0.0745     0.1156 1479.8794   -0.64   0.5194
## event1:extracted_run_number3    0.2402     0.1169 1479.8794    2.06   0.0400
##                                 
## (Intercept)                  ** 
## event1                       ** 
## extracted_run_number1        ***
## extracted_run_number2           
## extracted_run_number3        ***
## event1:extracted_run_number1    
## event1:extracted_run_number2    
## event1:extracted_run_number3 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 ext__1 ext__2 ext__3 e1:__1 e1:__2
## event1       0.000                                          
## extrctd_r_1 -0.001  0.000                                   
## extrctd_r_2 -0.001  0.000 -0.329                            
## extrctd_r_3  0.002  0.000 -0.338 -0.338                     
## evnt1:xt__1  0.000 -0.005  0.000  0.000  0.000              
## evnt1:xt__2  0.000 -0.005  0.000  0.000  0.000 -0.330       
## evnt1:xt__3  0.000  0.014  0.000  0.000  0.000 -0.337 -0.337
# checkmodel %>%
#   apa_print() %>%
#   apa_table()

plot(allEffects(checkmodel_bothdomains))

lsmeans(checkmodel_bothdomains, pairwise ~ event | extracted_run_number)$contrasts
## extracted_run_number = run1:
##  contrast    estimate    SE   df t.ratio p.value
##  unexp - exp    0.230 0.267 1480 0.860   0.3880 
## 
## extracted_run_number = run2:
##  contrast    estimate    SE   df t.ratio p.value
##  unexp - exp    0.286 0.267 1480 1.070   0.2840 
## 
## extracted_run_number = run3:
##  contrast    estimate    SE   df t.ratio p.value
##  unexp - exp    0.915 0.271 1480 3.380   0.0010 
## 
## extracted_run_number = run4:
##  contrast    estimate    SE   df t.ratio p.value
##  unexp - exp    0.308 0.267 1480 1.160   0.2480 
## 
## Degrees-of-freedom method: kenward-roger
model_habituation_results <- cbind(gen.m(checkmodel_bothdomains), gen.ci(checkmodel_bothdomains)[3:10,])
## Computing profile confidence intervals ...
model_habituation_results_byrun <- lsmeans(checkmodel_bothdomains, pairwise ~ event | extracted_run_number, pbkrtest.limit = 3048)$contrasts %>% as.data.frame()
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_psych-env_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_psych-env_perrun.Rds"))
univariate_data_both <- univariate_data %>%
    filter(event != "fam",
         # top_n_voxels == "100",
         domain == "both",
         str_length(extracted_run_number) == 4)

ROIs <- levels(as.factor(univariate_data_both$ROI_name_final))
runs <- c("run1", "run2", "run3", "run4", "all_runs", "runs12")
# top_n_voxels <- 100
modelsummaries_both_init <- data.frame()

5.2 Model tables

for (ROI in ROIs) {
  curROI <- ROI

  
  for (topN in top_n_voxels) {
    cur_top_voxels_N <- topN
    
    for (run in runs) {
      curRun <- run
      
      if (curRun != "all_runs" & curRun != "runs12") {
        ROIdata <- univariate_data_both %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN,
                 extracted_run_number == curRun)
        
      } else if (curRun == "all_runs") {
        ROIdata <- univariate_data_both %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN)
        
      } else if (curRun == "runs12") {
        ROIdata <- univariate_data_both %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN) %>%
          filter(extracted_run_number == "run1" |
                   extracted_run_number == "run2")
      }
      
      if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
        ROI_category <- "specific"
      }
      else {
        ROI_category <- "general"
      }
      
      model <- lmer(data = ROIdata,
                    formula = meanbeta ~ event + (1 | subjectID))
      
      focal_region <- ROIdata$focal_region[1]
      
      model_null <-
        update(model, formula = ~ . - event)  # Without interaction term
      
      BF_BIC_event <-
        data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
                   BF_interpretation = interpret_bf(exp((
                     BIC(model_null) - BIC(model)
                   ) / 2)))  # BICs to Bayes factor
      
      # BF_initial <- rbind(BF_BIC_event)
      
      rownames(BF_BIC_event) <- c("event1")
      
      BFs <- rbind(BF_intercept_row,
                   BF_BIC_event)
      
      effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
        select(d)
      effect_sizes <-
        rbind(effect_sizes_intercept_row, effect_sizes_initial)
      
      modelsummary <-
        cbind(
          cur_top_voxels_N,
          curRun,
          ROI_category,
          focal_region,
          curROI,
          rownames(summary(model)$coefficients),
          summary(model)$coefficients,
          confint.merMod(model)[3:4,],
          effect_sizes,
          BFs
        ) %>%
        as.data.frame()
      
      colnames(modelsummary) <-
        c(
          "top_n_voxels",
          "extracted_run_number",
          "domain",
          "focal_region",
          "ROI",
          "effect",
          "B",
          "SE",
          "df",
          "t",
          "p",
          "CI_95_lower",
          "CI_95_upper",
          "d",
          "BF",
          "BF_interpretation"
        )
      
      
      modelsummaries_both_init <-
        rbind(modelsummaries_both_init, modelsummary)
      cat(
        tab_model(
          model,
          show.std = TRUE,
          show.se = TRUE,
          show.stat = TRUE,
          auto.label = TRUE,
          title = paste0(
            "Exploratory Univariate ROI results for events from both domains, run = ",
            curRun,
            curROI
          )
        )$knitr,
        "  \n  \n"
      )
    }
  }
}
modelsummaries_both <- modelsummaries_both_init %>%
  mutate_at(c(7:15), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

modelsummaries_both$effect <- 
  as.factor(modelsummaries_both$effect)

levels(modelsummaries_both$effect) <-
  c("intercept", "event")

rownames(modelsummaries_both) <- NULL


# write_rds(modelsummaries_both, here("outputs/univariate_analysis_outputs/univar_results_psych-env_allregions_alltopN.Rds"))
modelsummaries_both <- read_rds( here("outputs/univariate_results/univar_results_psych-env_allregions_alltopN.Rds"))

5.3 Visualize

Effects by run and ROI size

modelsummaries_both$extracted_run_number <- factor(modelsummaries_both$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "runs12", "all_runs")) 

modelsummaries_both_focal <- modelsummaries_both %>% filter(focal_region == 1, effect == "event")

ggplot(modelsummaries_both_focal %>% filter(str_length(extracted_run_number) == 4), aes(extracted_run_number, d, colour = top_n_voxels)) +
  geom_point() +
  geom_line(aes(group = top_n_voxels)) +
  # geom_errorbar(aes(ymin = B - SE, ymax = B + SE, width = 0.1)) +
  # ylab("Bayes Factor (log)") +
  ylab("Cohen's D") +
  facet_wrap(~ effect + ROI, nrow = 1) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

5.3.1 Snsitivity to top voxel N

focal_modelsummaries <- modelsummaries_both %>%
  filter(focal_region == 1,
         str_length(extracted_run_number) == 4)

focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4"))

ggplot(
  focal_modelsummaries %>%
    filter(str_length(extracted_run_number) == 4,
           effect == "event"),
  aes(extracted_run_number, d, colour = top_n_voxels)
) +
  geom_point() +
  geom_line(aes(group = top_n_voxels)) +
  ylab("Cohen's D") +
  facet_wrap( ~ effect + factor(
    ROI,
    levels = c(
      "supramarginal_L",
      "supramarginal_R",
      "superiortemporal_L",
      "superiortemporal_R",
      "antParietal_bilateral",
      "precentral_inferiorfrontal_R",
      "V1_bilateral",
      "MT_bilateral"
    )
  ), nrow = 1) +
  # coord_cartesian(ylim = c(-.75, .75)) +
  geom_hline(yintercept = 0) +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  )) 

5.4 All regions, all runs

DT::datatable(modelsummaries_both %>% filter(extracted_run_number == "all_runs") %>% arrange(effect, focal_region), options = list(scrollX = TRUE))

5.5 All regions, runs 1-2

DT::datatable(modelsummaries_both %>% filter(extracted_run_number == "runs12") %>% arrange(effect, focal_region), options = list(scrollX = TRUE))

6 PART 4: Using VOE runs 3-4 to localize

Here we repeat the same analysis as confirmatory analyses above, except that we use VOE runs 3-4 to localize STS and SMG, and we search for voxels in the domain general ROI parcels that respond more to unexp than exp, to ask whether there are any voxels that show the VOE response in those regions, without considering the localizer data.

univariate_data_VOEextract <-
  readRDS(here("outputs", "univariate_data", "study2_univariate_data_VOEextract.Rds"))

univariate_summary_domain_VOEextract <-
  readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_VOEextract.Rds"))

univariate_summary_domain_runs12_VOEextract <-
  readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_runs12_VOEextract.Rds"))
focal_data_VOEextract <- univariate_data_VOEextract %>%
  filter(focal_region == 1)

focal_summary_domain_VOEextract <- univariate_summary_domain_VOEextract %>%
  filter(focal_region == 1)

focal_summary_domain_runs12_VOEextract <- univariate_summary_domain_runs12_VOEextract %>%
    filter(focal_region == 1)


domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")

regions <- levels(as.factor(focal_data_VOEextract$ROI_name_final))

6.1 Plots

First, here are plots of betas, run by run, to physics and psychology events (familiarization, expected, unexpected).

plot_univar_runbyrun("superiortemporal_L")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("superiortemporal_R")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("supramarginal_L")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("supramarginal_R")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("antParietal_bilateral")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("precentral_inferiorfrontal_R")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("V1_bilateral")
## Warning: Using alpha for a discrete variable is not advised.

plot_univar_runbyrun("MT_bilateral")
## Warning: Using alpha for a discrete variable is not advised.

DS_runs12_plots <- plot_univar_runs12("superiortemporal", 0, 4) |
plot_univar_runs12("supramarginal", 0, 4) 

DG_runs12_plots <- plot_univar_runs12("precentral_inferiorfrontal_R", 0, 4) | plot_univar_runs12("antParietal_bilateral", 0, 4) | plot_univar_runs12("V1_bilateral", 0, 11) | plot_univar_runs12("MT_bilateral", 0, 11)

DS_runs12_plots / DG_runs12_plots
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

focal_data_100_runs12_VOEextract <- focal_data_VOEextract %>%
  filter(
         event != "fam",
         top_n_voxels == "100",
         domain != "both",
         (extracted_run_number == "run1" | extracted_run_number == "run2"))

ROIs <- levels(as.factor(focal_data_100_runs12_VOEextract$ROI_name_final))

focal_data_100_runs12_VOEextract$event <- relevel(focal_data_100_runs12_VOEextract$event, ref = "exp")


focal_modelsummaries_100_VOEextract_init <- data.frame()

6.2 Model Tables

6.2.1 Psych-action and physics events

for (ROI in ROIs) {
  curROI <- ROI
  ROIdata <-
    focal_data_100_runs12_VOEextract %>% filter(ROI_name_final == curROI)
  # print("Current ROI: ", curROI)
  if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
    ROI_category <- "specific"
    model <- lmer(
      data = ROIdata,
      formula = meanbeta ~ event * domain + (1|subjectID)) # had to remove run random int (convergence))
      } else {
        ROI_category <- "general"
        model <- lmer(data = ROIdata,
                      formula = meanbeta ~ event * domain + (1|subjectID))
      }
  
  model_just_maineffects <- update(model, formula = ~ . -event:domain)  # Without interaction term
  model_just_domain <- update(model_just_maineffects, formula = ~ . -event)
  model_just_event <- update(model_just_maineffects, formula = ~ . -domain)
  
  BF_BIC_interaction <- data.frame(BF = exp((BIC(model_just_maineffects) - BIC(model))/2), BF_interpretation = interpret_bf(exp((BIC(model_just_maineffects) - BIC(model))/2)))  # BICs to Bayes factor
  
  BF_BIC_event <- data.frame(BF = exp((BIC(model_just_domain) - BIC(model_just_maineffects))/2),
                             BF_interpretation = interpret_bf(exp((BIC(model_just_domain) - BIC(model_just_maineffects))/2)))
  
  BF_BIC_domain <- data.frame(BF = exp((BIC(model_just_event) - BIC(model_just_maineffects))/2),
                              BF_interpretation = interpret_bf(exp((BIC(model_just_event) - BIC(model_just_maineffects))/2)))
  
  BF_initial <- rbind(BF_BIC_event,
                      BF_BIC_domain,
                      BF_BIC_interaction)
  
  rownames(BF_initial) <- rownames_no_intercept
  BFs <- rbind(BF_intercept_row,
               BF_initial)
    
  effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>% 
    select(d)
  effect_sizes <- rbind(effect_sizes_intercept_row, effect_sizes_initial)
  
  
      ROI_focal <- ROIdata$focal_region[1]

      modelsummary <-
        cbind(
          ROI_category,
          ROI_focal,
          curROI,
          rownames(summary(model)$coefficients),
          summary(model)$coefficients,
          confint.merMod(model)[3:6, ],
          effect_sizes,
          BFs
        ) %>%
        as.data.frame()
      
      colnames(modelsummary) <-
        c(
          "domain",
          "focal_region",
          "ROI",
          "effect",
          "B",
          "SE",
          "df",
          "t",
          "p",
          "CI_95_lower",
          "CI_95_upper",
          "d",
          "BF",
          "BF_interpretation"
        )
 
  focal_modelsummaries_100_VOEextract_init <-
    rbind(focal_modelsummaries_100_VOEextract_init, modelsummary)
  # model_parameters(model) %>% print_md()
  cat(
    tab_model(
      model,
      show.std = TRUE,
      show.se = TRUE,
      show.stat = TRUE,
      auto.label = TRUE,
      title = paste0("Explratory Univariate ROI results (localized using VOE runs 3-4) in: ", curROI)
    )$knitr,
    "  \n  \n"
  )
  # focal_modeltables_100 <- append(focal_modeltables_100, cur_tab)
  
  }
focal_modelsummaries_100_VOEextract <- focal_modelsummaries_100_VOEextract_init %>%
  mutate_at(c(5:13), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

focal_modelsummaries_100_VOEextract$effect <-
  as.factor(focal_modelsummaries_100_VOEextract$effect)

levels(focal_modelsummaries_100_VOEextract$effect) <-
  c("intercept", "domain", "event", "event:domain")

rownames(focal_modelsummaries_100_VOEextract) <- NULL


# write_rds(focal_modelsummaries_100_VOEextract, here("outputs/univariate_results/univar_results_top100_VOE-extract_runs12.Rds"))
focal_modelsummaries_100_VOEextract <- read_rds(here("outputs/univariate_results/univar_results_top100_VOE-extract_runs12.Rds"))

DT::datatable(dplyr::arrange(focal_modelsummaries_100_VOEextract, effect), options = list(scrollX = TRUE))

6.2.2 Psych-env events

focal_data_100_runs1234_VOEextract_bothdomains <-
  focal_data_VOEextract %>%
  filter(
    event != "fam",
    top_n_voxels == "100",
    domain == "both",
    str_length(extracted_run_number) == 4
  )

ROIs <-
  levels(as.factor(
    focal_data_100_runs1234_VOEextract_bothdomains$ROI_name_final
  ))

focal_data_100_runs1234_VOEextract_bothdomains$event <-
  relevel(focal_data_100_runs1234_VOEextract_bothdomains$event, ref = "unexp")


focal_modelsummaries_100_VOEextract_bothdomains_init <- data.frame()
for (ROI in ROIs) {
  curROI <- ROI
  cur_top_voxels_N <- 100  
  curRun <- "runs12"

  ROIdata <- focal_data_100_runs1234_VOEextract_bothdomains %>%
    filter(ROI_name_final == curROI) %>%
    filter(extracted_run_number == "run1" |
             extracted_run_number == "run2")
  
  if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
    ROI_category <- "specific"
  }
  else {
    ROI_category <- "general"
  }
  
  model <- lmer(data = ROIdata,
                formula = meanbeta ~ event + (1 | subjectID))
  
  focal_region <- ROIdata$focal_region[1]
  
  model_null <-
    update(model, formula = ~ . - event)  # Without interaction term
  
  BF_BIC_event <-
    data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
               BF_interpretation = interpret_bf(exp((
                 BIC(model_null) - BIC(model)
               ) / 2)))  # BICs to Bayes factor
  
  # BF_initial <- rbind(BF_BIC_event)
  
  rownames(BF_BIC_event) <- c("event1")
  
  BFs <- rbind(BF_intercept_row,
               BF_BIC_event)
  
  effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
    select(d)
  effect_sizes <-
    rbind(effect_sizes_intercept_row, effect_sizes_initial)
  
  modelsummary <-
    cbind(
      cur_top_voxels_N,
      curRun,
      ROI_category,
      focal_region,
      curROI,
      rownames(summary(model)$coefficients),
      summary(model)$coefficients,
      confint.merMod(model)[3:4, ],
      effect_sizes,
      BFs
    ) %>%
    as.data.frame()
  
  colnames(modelsummary) <-
    c(
      "top_n_voxels",
      "extracted_run_number",
      "domain",
      "focal_region",
      "ROI",
      "effect",
      "B",
      "SE",
      "df",
      "t",
      "p",
      "CI_95_lower",
      "CI_95_upper",
      "d",
      "BF",
      "BF_interpretation"
    )
      
  
  
  focal_modelsummaries_100_VOEextract_bothdomains_init <-
    rbind(focal_modelsummaries_100_VOEextract_bothdomains_init,
          modelsummary)
  cat(
    tab_model(
      model,
      show.std = TRUE,
      show.se = TRUE,
      show.stat = TRUE,
      auto.label = TRUE,
      title = paste0(
        "Exploratory Univariate ROI results for events from both domains, run = ",
        curRun,
        curROI
      )
    )$knitr,
    "  \n  \n"
  )
}
focal_modelsummaries_100_VOEextract_bothdomains <- focal_modelsummaries_100_VOEextract_bothdomains_init %>%
  mutate_at(c(7:15), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

focal_modelsummaries_100_VOEextract_bothdomains$effect <- 
  as.factor(focal_modelsummaries_100_VOEextract_bothdomains$effect)

levels(focal_modelsummaries_100_VOEextract_bothdomains$effect) <-
  c("intercept", "event")

rownames(focal_modelsummaries_100_VOEextract_bothdomains) <- NULL

# write_rds(focal_modelsummaries_100_VOEextract_bothdomains, here("outputs/univariate_results/univar_results_psych-env_top100_VOE-extract_runs12.Rds"))
focal_modelsummaries_100_VOEextract_bothdomains <- read_rds( here("outputs/univariate_results/univar_results_psych-env_top100_VOE-extract_runs12.Rds"))
DT::datatable(focal_modelsummaries_100_VOEextract_bothdomains, options = list(scrollX = TRUE))

7 PART 5: Single betas

single_betas0 <- readRDS(here("outputs/univariate_data/study2_univariate_data_singlebetas.Rds")) %>%
  select(-ROI_category)

regioninfo_category <- regioninfo %>%
  select(ROI_name_final, ROI_category)

single_betas <- left_join(single_betas0, regioninfo_category)
## Joining, by = "ROI_name_final"
single_betas$ROI_name_final <- as.factor(single_betas$ROI_name_final) 

single_betas$event <- as.factor(single_betas$event)
single_betas$event <- relevel(single_betas$event, ref = "fam")
single_betas$task <- factor(single_betas$task, levels = c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint"))
single_betas_summary_domain <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_domain_singlebetas.Rds"))

single_betas_summary_task <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_task_singlebetas.Rds"))

single_betas_summary_scenario <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_scenario_singlebetas.Rds"))

single_betas_summary_scenario$event <- relevel(single_betas_summary_scenario$event, ref = "fam")
single_betas_focal <- single_betas %>%
  filter(focal_region == 1)

7.1 Plots

# solidity - 25a88b
# permanence - 19626a
# goal - f24b00
# efficiency - deae1e
# agent-solidity - 6600cd 
# missing-barrier - fba0d0

  
ggplot(
  single_betas_focal %>%
    filter(
      str_length(extracted_run_number) == 4,
      event != "fam"
    ),
  aes(task, meanbeta, alpha = event, fill = task)
) +
  geom_bar(
    position = "dodge",
    stat = "summary",
    fun.data = "mean_se",
    colour = "black"
  ) +
  geom_errorbar(
    position = "dodge",
    stat = "summary",
    fun.data = "mean_se",
    colour = "black"
  ) +
  scale_fill_manual(values = c(
    "#25a88b",
    "#19626a",
    "#f24b00",
    "#deae1e",
    "#6600cd",
    "#fba0d0"
  )) + # ,
  # geom_point() +
  theme(axis.text.x = element_blank()) +
  # scale_alpha_continuous(range = c(0.5, 1)) +
  facet_wrap( ~ ROI_name_final, scales = "free", nrow = 2) +
  ggtitle("Runs 1-4, Focal regions")
## Warning: Using alpha for a discrete variable is not advised.

ggplot(
  single_betas_focal %>%
    filter(ROI_name_final == "supramarginal_L"),
  aes(
    vis_statistics_video_group,
    meanbeta,
    alpha = event,
    fill = task
  )
) +
  geom_bar(
    position = "dodge",
    stat = "summary",
    fun.data = "mean_se",
    colour = "black"
  ) +
  geom_errorbar(
    position = "dodge",
    stat = "summary",
    fun.data = "mean_se",
    colour = "black"
  ) +
  scale_fill_manual(values = c(
    "#25a88b",
    "#19626a",
    "#f24b00",
    "#deae1e",
    "#6600cd",
    "#fba0d0"
  )) +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  )) +
  ggtitle("Region: supramarginal_L")
## Warning: Using alpha for a discrete variable is not advised.

Plotting each event type in each task against run number, per ROI.

ggplot(single_betas_focal, aes(as.integer(extracted_run_number), meanbeta, colour = event)) +
  geom_smooth(method = "glm") +
  facet_wrap( ~ ROI_name_final + task, nrow = 3, scales = "free") +
  scale_colour_brewer(palette = "Set2")
## `geom_smooth()` using formula = 'y ~ x'

Visualizing visual features against responses.

ggplot(single_betas_focal,
       aes(normalized_iqr_contrast, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(single_betas_focal,
       aes(normalized_mean_luminance, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(single_betas_focal,
       aes(normalized_mean_motion, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(single_betas_focal,
       aes(normalized_mean_hsf, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(single_betas_focal,
       aes(normalized_mean_lsf, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(single_betas_focal,
       aes(normalized_mean_rect, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(single_betas_focal,
       aes(normalized_mean_curv, meanbeta)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

7.2 Modeling

single_betas_focal$event <- relevel(single_betas_focal$event, ref = "unexp")
single_betas_focal_psychphys_runs12 <- single_betas_focal %>%
  filter(domain != "psychology-environment",
         event != "fam") %>%
  filter(extracted_run_number == "run1" | extracted_run_number == "run2")

ROIs <- unique(single_betas_focal_psychphys_runs12$ROI_name_final) %>% as.vector()

single_betas_focal_psychphys_runs12_summaries <- data.frame()
rownames_no_intercept <- c("event1", "domain1", "event:domain")

7.3 Model Tables

for (ROI in ROIs) {
  curROI <- ROI
  ROIdata <-
    single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == curROI)
  # print("Current ROI: ", curROI)
  if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
    ROI_category <- "specific"
    model <- lmer(
      data = ROIdata,
      formula = meanbeta ~ event * domain + (1|subjectID) + (1|vis_statistics_video_group) 
    ) # had to remove run random int (convergence))
  } else {
    ROI_category <- "general"
    model <- lmer(
      data = ROIdata,
      formula = meanbeta ~ event * domain + (1|subjectID) + (1|vis_statistics_video_group)
    )
  }
  
  model_just_maineffects <-
    update(model, formula = ~ . - event:domain)  # Without interaction term
  model_just_domain <-
    update(model_just_maineffects, formula = ~ . - event)
  model_just_event <-
    update(model_just_maineffects, formula = ~ . - domain)
  
  BF_BIC_interaction <-
    data.frame(BF = exp((
      BIC(model_just_maineffects) - BIC(model)
    ) / 2), BF_interpretation = interpret_bf(exp((
      BIC(model_just_maineffects) - BIC(model)
    ) / 2)))  # BICs to Bayes factor
  
  BF_BIC_event <-
    data.frame(BF = exp((
      BIC(model_just_domain) - BIC(model_just_maineffects)
    ) / 2),
    BF_interpretation = interpret_bf(exp((
      BIC(model_just_domain) - BIC(model_just_maineffects)
    ) / 2)))
  
  BF_BIC_domain <-
    data.frame(BF = exp((
      BIC(model_just_event) - BIC(model_just_maineffects)
    ) / 2),
    BF_interpretation = interpret_bf(exp((
      BIC(model_just_event) - BIC(model_just_maineffects)
    ) / 2)))
  
  BF_initial <- rbind(BF_BIC_event,
                      BF_BIC_domain,
                      BF_BIC_interaction)
  
  rownames(BF_initial) <- rownames_no_intercept
  BFs <- rbind(BF_intercept_row,
               BF_initial)
  
  effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
    select(d)
  effect_sizes <-
    rbind(effect_sizes_intercept_row, effect_sizes_initial)
  ROI_focal <- ROIdata$focal_region[1]
  cur_top_voxels_N <- 100
  curRun <- "runs12"
  
  modelsummary <-
    cbind(
      cur_top_voxels_N,
      curRun,
      ROI_category,
      ROI_focal,
      curROI,
      rownames(summary(model)$coefficients),
      summary(model)$coefficients,
      confint.merMod(model)[3:6,],
      effect_sizes,
      BFs
    ) %>%
    as.data.frame()
  
  colnames(modelsummary) <-
    c(
      "top_n_voxels",
      "extracted_run_number",
      "domain",
      "focal_region",
      "ROI",
      "effect",
      "B",
      "SE",
      "df",
      "t",
      "p",
      "CI_95_lower",
      "CI_95_upper",
      "d",
      "BF",
      "BF_interpretation"
    )
  
  single_betas_focal_psychphys_runs12_summaries <-
    rbind(single_betas_focal_psychphys_runs12_summaries,
          modelsummary)
  # model_parameters(model) %>% print_md()
  cat(
    tab_model(
      model,
      show.std = TRUE,
      show.se = TRUE,
      show.stat = TRUE,
      auto.label = TRUE,
      title = paste0("Single beta Univariate ROI results (runs 1-2) in: ", curROI)
    )$knitr,
    "  \n  \n"
  )
  # focal_modeltables_100 <- append(focal_modeltables_100, cur_tab)
  
}
single_betas_focal_psychphys_runs12_summaries <- single_betas_focal_psychphys_runs12_summaries %>%
  mutate_at(c(7:15), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

single_betas_focal_psychphys_runs12_summaries$effect <-
  as.factor(single_betas_focal_psychphys_runs12_summaries$effect)


levels(single_betas_focal_psychphys_runs12_summaries$effect) <-
  c("intercept", "domain", "event", "event:domain")

rownames(single_betas_focal_psychphys_runs12_summaries) <- NULL


# write_rds(single_betas_focal_psychphys_runs12_summaries, here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys.Rds"))
single_betas_focal_psychphys_runs12_summaries <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys.Rds"))

DT::datatable(dplyr::arrange(single_betas_focal_psychphys_runs12_summaries, effect), options = list(scrollX = TRUE))

7.4 PART 5a: Visual statistics (psych and physics)

single_betas_focal_visualstats_summary_init <- data.frame()


for (ROI in ROIs) {
  curROI <- ROI
  ROIdata <-
    single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == curROI)
  # print("Current ROI: ", curROI)
  if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
    ROI_category <- "specific"
      } else {
        ROI_category <- "general"
      }
    model <- lmer(data = ROIdata, formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
      effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
        select(d)
      effect_sizes <-
        rbind(effect_sizes_intercept_row, effect_sizes_initial)
      
  ROI_focal <- ROIdata$focal_region[1]
  cur_top_voxels_N <- 100
  curRun <- "runs12"
  
  modelsummary <-
    cbind(
      cur_top_voxels_N,
      curRun,
      ROI_category,
      ROI_focal,
      curROI,
      rownames(summary(model)$coefficients),
      summary(model)$coefficients,
      confint.merMod(model)[3:12,],
      effect_sizes
      # BFs
    ) %>%
    as.data.frame()
  
  colnames(modelsummary) <-
    c(
      "top_n_voxels",
      "extracted_run_number",
      "domain",
      "focal_region",
      "ROI",
      "effect",
      "B",
      "SE",
      "df",
      "t",
      "p",
      "CI_95_lower",
      "CI_95_upper",
      "d"
    )
  
    single_betas_focal_visualstats_summary_init <-
    rbind(single_betas_focal_visualstats_summary_init,
          modelsummary)
    
  # model_parameters(model) %>% print_md()
  # cat(
  #   tab_model(
  #     model,
  #     show.std = TRUE,
  #     show.se = TRUE,
  #     show.stat = TRUE,
  #     auto.label = TRUE,
  #     title = paste0("Exploratory Univariate ROI results in: ", curROI)
  #   )$knitr,
  #   "  \n  \n"
  # )
  }
single_betas_focal_visualstats_summary <- single_betas_focal_visualstats_summary_init %>%
  mutate_at(c(7:14), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

single_betas_focal_visualstats_summary$effect <-
  as.factor(single_betas_focal_visualstats_summary$effect)


levels(single_betas_focal_visualstats_summary$effect)[1:4] <-
  c("intercept", "domain", "event", "event:domain")

rownames(single_betas_focal_visualstats_summary) <- NULL


# write_rds(single_betas_focal_visualstats_summary, here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys_visstats.Rds"))
single_betas_focal_visualstats_summary <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys_visstats.Rds"))

DT::datatable(dplyr::arrange(single_betas_focal_visualstats_summary, effect), options = list(scrollX = TRUE))
V1_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "V1_bilateral"),
                             formula = meanbeta ~ event * domain +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf  + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID)) 

summary(V1_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | extracted_run_number) + (1 |      subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "V1_bilateral")
## 
## REML criterion at convergence: 3916
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.065 -0.581 -0.022  0.648  4.223 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 12.5313  3.540   
##  extracted_run_number (Intercept)  0.0802  0.283   
##  Residual                          8.1251  2.850   
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 9.0972     0.6683  25.6178   13.61  3.1e-13 ***
## event1                      0.0655     0.1033 726.0307    0.63   0.5259    
## domain1                     0.4682     0.1767 726.0014    2.65   0.0082 ** 
## normalized_iqr_contrast    -0.0388     0.1052 726.0004   -0.37   0.7119    
## normalized_mean_luminance  -1.2349     0.2240 726.0008   -5.51  4.9e-08 ***
## normalized_mean_motion      0.7183     0.1649 726.0038    4.36  1.5e-05 ***
## normalized_mean_hsf         0.3798     0.2439 726.0013    1.56   0.1198    
## normalized_mean_rect       -0.3466     0.2017 726.0002   -1.72   0.0861 .  
## normalized_mean_curv        0.1841     0.1807 726.0021    1.02   0.3086    
## event1:domain1              0.0613     0.1030 726.0002    0.60   0.5520    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.003                                                            
## domain1     -0.003 -0.055                                                     
## nrmlzd_qr_c -0.062 -0.036  0.084                                              
## nrmlzd_mn_l  0.048 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.040  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.005 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.016 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.042 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.002 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
check_collinearity(V1_model_singlebetas)
MT_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "MT_bilateral"),
                             formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID)) 

summary(MT_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | extracted_run_number) + (1 |      subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "MT_bilateral")
## 
## REML criterion at convergence: 3098
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.229 -0.591  0.009  0.605  3.287 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 6.351    2.520   
##  extracted_run_number (Intercept) 0.016    0.126   
##  Residual                         2.716    1.648   
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                            Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                 5.42973    0.45984  30.85320   11.81  5.6e-13 ***
## event1                      0.04583    0.05970 726.04405    0.77  0.44295    
## domain1                     0.18927    0.10213 726.00118    1.85  0.06427 .  
## normalized_iqr_contrast    -0.53523    0.06080 725.99966   -8.80  < 2e-16 ***
## normalized_mean_luminance  -0.25298    0.12949 726.00030   -1.95  0.05112 .  
## normalized_mean_motion      0.35157    0.09535 726.00465    3.69  0.00024 ***
## normalized_mean_hsf        -1.38930    0.14100 726.00108   -9.85  < 2e-16 ***
## normalized_mean_rect        0.23006    0.11660 725.99933    1.97  0.04887 *  
## normalized_mean_curv       -0.59928    0.10449 726.00211   -5.74  1.4e-08 ***
## event1:domain1              0.00941    0.05954 725.99932    0.16  0.87444    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.002                                                            
## domain1     -0.003 -0.055                                                     
## nrmlzd_qr_c -0.052 -0.036  0.084                                              
## nrmlzd_mn_l  0.040 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.034  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.004 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.013 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.035 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.002 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
check_collinearity(MT_model_singlebetas)
APC_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "antParietal_bilateral"),
                             formula = meanbeta ~ event * domain +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf  + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID)) 

summary(APC_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | extracted_run_number) + (1 |      subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "antParietal_bilateral")
## 
## REML criterion at convergence: 3635
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.046 -0.535  0.015  0.520  4.165 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 4.3797   2.093   
##  extracted_run_number (Intercept) 0.0278   0.167   
##  Residual                         5.7731   2.403   
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 3.7034     0.4017  23.0478    9.22  3.4e-09 ***
## event1                      0.2533     0.0870 726.0513    2.91  0.00372 ** 
## domain1                     0.1349     0.1489 726.0023    0.91  0.36517    
## normalized_iqr_contrast    -0.2247     0.0887 726.0005   -2.54  0.01145 *  
## normalized_mean_luminance  -0.3695     0.1888 726.0013   -1.96  0.05073 .  
## normalized_mean_motion      0.2846     0.1390 726.0063    2.05  0.04095 *  
## normalized_mean_hsf        -0.6989     0.2056 726.0022   -3.40  0.00071 ***
## normalized_mean_rect       -0.0208     0.1700 726.0002   -0.12  0.90257    
## normalized_mean_curv       -0.4586     0.1523 726.0034   -3.01  0.00270 ** 
## event1:domain1             -0.0903     0.0868 726.0002   -1.04  0.29868    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.004                                                            
## domain1     -0.005 -0.055                                                     
## nrmlzd_qr_c -0.087 -0.036  0.084                                              
## nrmlzd_mn_l  0.067 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.056  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.007 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.022 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.058 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.003 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
check_collinearity(APC_model_singlebetas)
RFC_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "precentral_inferiorfrontal_R"),
                             formula = meanbeta ~ event * domain +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID)) 

summary(RFC_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "precentral_inferiorfrontal_R")
## 
## REML criterion at convergence: 3629
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.369 -0.594  0.030  0.630  3.362 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 3.19     1.79    
##  Residual              5.81     2.41    
## Number of obs: 768, groups:  subjectID, 32
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 3.8432     0.3324  32.8136   11.56  4.1e-13 ***
## event1                      0.2549     0.0873 727.0000    2.92   0.0036 ** 
## domain1                     0.2856     0.1493 727.0000    1.91   0.0562 .  
## normalized_iqr_contrast    -0.2156     0.0889 727.0000   -2.43   0.0155 *  
## normalized_mean_luminance  -0.2577     0.1893 727.0000   -1.36   0.1738    
## normalized_mean_motion      0.1398     0.1394 727.0000    1.00   0.3161    
## normalized_mean_hsf        -0.6089     0.2062 727.0000   -2.95   0.0032 ** 
## normalized_mean_rect        0.1171     0.1705 727.0000    0.69   0.4926    
## normalized_mean_curv       -0.3032     0.1528 727.0000   -1.98   0.0475 *  
## event1:domain1             -0.0143     0.0871 727.0000   -0.16   0.8696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.005                                                            
## domain1     -0.006 -0.055                                                     
## nrmlzd_qr_c -0.106 -0.036  0.084                                              
## nrmlzd_mn_l  0.081 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.068  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.009 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.027 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.071 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.003 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
check_collinearity(RFC_model_singlebetas)
RSTS_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "superiortemporal_R"),
                             formula = meanbeta ~ event * domain +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf  + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID)) 

summary(RSTS_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | extracted_run_number) + (1 |      subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "superiortemporal_R")
## 
## REML criterion at convergence: 3402
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.790 -0.573  0.014  0.589  3.773 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 2.50595  1.5830  
##  extracted_run_number (Intercept) 0.00173  0.0416  
##  Residual                         4.29415  2.0722  
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 2.0826     0.2951  28.1989    7.06  1.1e-07 ***
## event1                      0.1045     0.0751 726.1236    1.39  0.16439    
## domain1                    -0.3892     0.1284 726.0052   -3.03  0.00253 ** 
## normalized_iqr_contrast    -0.1251     0.0765 726.0009   -1.64  0.10224    
## normalized_mean_luminance  -0.0367     0.1628 726.0027   -0.23  0.82180    
## normalized_mean_motion     -0.0225     0.1199 726.0149   -0.19  0.85101    
## normalized_mean_hsf        -0.6320     0.1773 726.0049   -3.56  0.00039 ***
## normalized_mean_rect        0.1491     0.1466 726.0000    1.02  0.30962    
## normalized_mean_curv       -0.0527     0.1314 726.0078   -0.40  0.68854    
## event1:domain1              0.0298     0.0749 725.9999    0.40  0.69066    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.005                                                            
## domain1     -0.005 -0.055                                                     
## nrmlzd_qr_c -0.103 -0.036  0.084                                              
## nrmlzd_mn_l  0.078 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.066  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.009 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.026 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.069 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.003 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
check_collinearity(RSTS_model_singlebetas)
LSTS_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "superiortemporal_L"),
                             formula = meanbeta ~ event * domain + event_n_within_run +domain * event +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf  + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID)) 

summary(LSTS_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + event_n_within_run + domain * event +  
##     normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion +  
##     normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv +  
##     (1 | extracted_run_number) + (1 | subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "superiortemporal_L")
## 
## REML criterion at convergence: 3719
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.839 -0.511  0.058  0.546  3.911 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 4.386    2.094   
##  extracted_run_number (Intercept) 0.111    0.333   
##  Residual                         6.414    2.533   
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)                 1.1481     0.4877  11.9224    2.35    0.037 * 
## event1                      0.0625     0.0919 725.0360    0.68    0.496   
## domain1                    -0.3916     0.1570 724.9993   -2.50    0.013 * 
## event_n_within_run         -0.0305     0.0102 732.2838   -2.98    0.003 **
## normalized_iqr_contrast    -0.0894     0.0934 724.9989   -0.96    0.339   
## normalized_mean_luminance   0.1916     0.1990 725.0016    0.96    0.336   
## normalized_mean_motion     -0.0932     0.1465 725.0002   -0.64    0.525   
## normalized_mean_hsf        -0.3465     0.2167 725.0001   -1.60    0.110   
## normalized_mean_rect       -0.2323     0.1792 725.0003   -1.30    0.195   
## normalized_mean_curv        0.0402     0.1607 725.0061    0.25    0.803   
## event1:domain1             -0.0997     0.0916 725.0082   -1.09    0.276   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 evn___ nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m
## event1      -0.015                                                       
## domain1      0.000 -0.055                                                
## evnt_n_wth_ -0.376  0.050 -0.010                                         
## nrmlzd_qr_c -0.072 -0.037  0.084 -0.011                                  
## nrmlzd_mn_l  0.066 -0.027  0.212 -0.021 -0.032                           
## nrmlzd_mn_m  0.047  0.086 -0.585  0.004 -0.381    -0.212                 
## nrmlzd_mn_h -0.001 -0.062  0.714 -0.014  0.399     0.344      -0.605     
## nrmlzd_mn_r  0.013 -0.009  0.324  0.018 -0.017    -0.432      -0.089     
## nrmlzd_mn_c -0.039 -0.056  0.442 -0.031  0.365     0.498      -0.615     
## event1:dmn1  0.011 -0.006  0.030 -0.037  0.017     0.002      -0.049     
##             nrmlzd_mn_h nrmlzd_mn_r nrmlzd_mn_c
## event1                                         
## domain1                                        
## evnt_n_wth_                                    
## nrmlzd_qr_c                                    
## nrmlzd_mn_l                                    
## nrmlzd_mn_m                                    
## nrmlzd_mn_h                                    
## nrmlzd_mn_r  0.296                             
## nrmlzd_mn_c  0.497      -0.426                 
## event1:dmn1  0.030       0.011       0.024
check_collinearity(LSTS_model_singlebetas)
RSMG_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "supramarginal_R"),
                             formula = meanbeta ~ event * domain +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv +  (1|extracted_run_number) + (1|subjectID)) 

summary(RSMG_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | extracted_run_number) + (1 |      subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "supramarginal_R")
## 
## REML criterion at convergence: 3449
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.579 -0.562  0.005  0.565  6.309 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 7.3025   2.702   
##  extracted_run_number (Intercept) 0.0787   0.281   
##  Residual                         4.3692   2.090   
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 2.5859     0.5250  20.0189    4.93 0.000081 ***
## event1                      0.1017     0.0757 726.0185    1.34    0.179    
## domain1                     0.2934     0.1295 726.0007    2.27    0.024 *  
## normalized_iqr_contrast    -0.1754     0.0771 726.0001   -2.27    0.023 *  
## normalized_mean_luminance  -0.4043     0.1642 726.0003   -2.46    0.014 *  
## normalized_mean_motion      0.1410     0.1209 726.0021    1.17    0.244    
## normalized_mean_hsf        -0.4244     0.1788 726.0007   -2.37    0.018 *  
## normalized_mean_rect        0.1570     0.1479 725.9999    1.06    0.289    
## normalized_mean_curv       -0.2804     0.1325 726.0011   -2.12    0.035 *  
## event1:domain1              0.1810     0.0755 725.9999    2.40    0.017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.003                                                            
## domain1     -0.003 -0.055                                                     
## nrmlzd_qr_c -0.058 -0.036  0.084                                              
## nrmlzd_mn_l  0.044 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.037  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.005 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.015 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.039 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.002 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
plot(allEffects(RSMG_model_singlebetas))

lsmeans(RSMG_model_singlebetas, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.565 0.213 726  2.649  0.0080 
## 
## domain = psychology-action:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp   -0.158 0.214 726 -0.739  0.4600 
## 
## Degrees-of-freedom method: kenward-roger
check_collinearity(RSMG_model_singlebetas)
LSMG_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
                             formula = meanbeta ~ event * domain +  normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf  + normalized_mean_rect + normalized_mean_curv +  (1|extracted_run_number) + (1|subjectID)) 

summary(LSMG_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +  
##     normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +  
##     normalized_mean_curv + (1 | extracted_run_number) + (1 |      subjectID)
##    Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==  
##     "supramarginal_L")
## 
## REML criterion at convergence: 3379
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.060 -0.600  0.012  0.579  5.640 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  subjectID            (Intercept) 4.11543  2.0287  
##  extracted_run_number (Intercept) 0.00379  0.0616  
##  Residual                         4.07805  2.0194  
## Number of obs: 768, groups:  subjectID, 32; extracted_run_number, 2
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 2.4037     0.3715  30.2024    6.47  3.7e-07 ***
## event1                      0.0665     0.0731 726.1058    0.91   0.3636    
## domain1                     0.2683     0.1252 726.0045    2.14   0.0324 *  
## normalized_iqr_contrast    -0.0619     0.0745 726.0008   -0.83   0.4064    
## normalized_mean_luminance  -0.4990     0.1587 726.0024   -3.14   0.0017 ** 
## normalized_mean_motion     -0.0931     0.1168 726.0128   -0.80   0.4256    
## normalized_mean_hsf        -0.0259     0.1728 726.0042   -0.15   0.8808    
## normalized_mean_rect        0.2933     0.1429 726.0001    2.05   0.0405 *  
## normalized_mean_curv       -0.0639     0.1280 726.0067   -0.50   0.6176    
## event1:domain1              0.2042     0.0730 726.0000    2.80   0.0053 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1       0.004                                                            
## domain1     -0.004 -0.055                                                     
## nrmlzd_qr_c -0.079 -0.036  0.084                                              
## nrmlzd_mn_l  0.061 -0.026  0.212 -0.033                                       
## nrmlzd_mn_m  0.051  0.086 -0.585 -0.381    -0.212                             
## nrmlzd_mn_h -0.007 -0.062  0.714  0.399     0.344      -0.605                 
## nrmlzd_mn_r  0.020 -0.010  0.325 -0.017    -0.432      -0.089       0.296     
## nrmlzd_mn_c -0.053 -0.054  0.442  0.365     0.497      -0.615       0.497     
## event1:dmn1 -0.003 -0.004  0.029  0.017     0.002      -0.049       0.029     
##             nrmlzd_mn_r nrmlzd_mn_c
## event1                             
## domain1                            
## nrmlzd_qr_c                        
## nrmlzd_mn_l                        
## nrmlzd_mn_m                        
## nrmlzd_mn_h                        
## nrmlzd_mn_r                        
## nrmlzd_mn_c -0.425                 
## event1:dmn1  0.012       0.023
plot(allEffects(LSMG_model_singlebetas))

lsmeans(LSMG_model_singlebetas, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.541 0.206 726  2.626  0.0088 
## 
## domain = psychology-action:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp   -0.275 0.207 726 -1.330  0.1838 
## 
## Degrees-of-freedom method: kenward-roger
check_collinearity(LSMG_model_singlebetas)

7.5 PART 5b: Per task event analysis

single_betas_pertask_data <- single_betas %>%
  filter(event != "fam") %>%
  filter(str_length(extracted_run_number) == 4)

tasks <- unique(single_betas_pertask_data$task)
ROIs <- unique(single_betas_pertask_data$ROI_name_final)
single_betas_pertask_data$event <- relevel(single_betas_pertask_data$event, ref = "unexp")

runs <- c("all_runs", "runs12")
single_betas_pertask_results <- data.frame()


for (ROI in ROIs) {
  curROI <- ROI
  for (task in tasks) {
    cur_task <- task
    
    for (run in runs) {
      if (run == "all_runs") {
        ROIdata <-
          single_betas_pertask_data %>% filter(ROI_name_final == curROI,
                                               task == cur_task)
      }
      
      else if (run == "runs12") {
        ROIdata <-
          single_betas_pertask_data %>%
          filter(ROI_name_final == curROI,
                 task == cur_task) %>%
          filter(extracted_run_number == "run1" |
                   extracted_run_number == "run2")
        
      }
      # print("Current ROI: ", curROI)
      
      
      model <- lmer(data = ROIdata,
                    formula = meanbeta ~ event + (1 | subjectID)) # had to remove run random int (convergence))
      
      focal_region <- ROIdata$focal_region[1]
      task_domain <- as.character(ROIdata$domain[1])
      ROI_category <- ROIdata$ROI_category[1]
      
      
      model_null <-
        update(model, formula = ~ . - event)  # Without interaction term
      
      BF_BIC_event <-
        data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
                   BF_interpretation = interpret_bf(exp((
                     BIC(model_null) - BIC(model)
                   ) / 2)))  # BICs to Bayes factor
      
      # BF_initial <- rbind(BF_BIC_event)
      
      rownames(BF_BIC_event) <- c("event1")
      
      BFs <- rbind(BF_intercept_row,
                   BF_BIC_event)
      
      effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
        select(d)
      effect_sizes <-
        rbind(effect_sizes_intercept_row, effect_sizes_initial)
      cur_top_voxels_N <- 100
      curRun <- run
      modelsummary <-
        cbind(
          cur_top_voxels_N,
          curRun,
          ROI_category,
          focal_region,
          curROI,
          cur_task,
          rownames(summary(model)$coefficients),
          summary(model)$coefficients,
          confint.merMod(model)[3:4, ],
          effect_sizes,
          BFs
        ) %>%
        as.data.frame()
      
      colnames(modelsummary) <-
        c(
          "top_n_voxels",
          "extracted_run_number",
          "domain",
          "focal_region",
          "ROI",
          "task",
          "effect",
          "B",
          "SE",
          "df",
          "t",
          "p",
          "CI_95_lower",
          "CI_95_upper",
          "d",
          "BF",
          "BF_interpretation"
        )
      
      
      single_betas_pertask_results <-
        rbind(single_betas_pertask_results, modelsummary)
      # model_parameters(model) %>% print_md()
      cat(
        tab_model(
          model,
          show.std = TRUE,
          show.se = TRUE,
          show.stat = TRUE,
          auto.label = TRUE,
          title = paste0(
            "Exploratory Univariate ROI results for the task ",
            cur_task,
            "in the region: ",
            curROI,
            "in runs ",
            run
          )
        )$knitr,
        "  \n  \n"
      )
    }
  }
}
single_betas_pertask_results <- single_betas_pertask_results %>%
  mutate_at(c(8:16), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

single_betas_pertask_results$effect <-
  as.factor(single_betas_pertask_results$effect)


levels(single_betas_pertask_results$effect) <-
  c("intercept", "event")

rownames(single_betas_pertask_results) <- NULL


# write_rds(single_betas_pertask_results, here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
single_betas_pertask_results <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
DT::datatable(dplyr::arrange(single_betas_pertask_results, extracted_run_number, effect), options = list(scrollX = TRUE))

8 Part 6: Overlap of top 100 voxels

overlap_data_prelim <- read_csv(here("input_data/ROI_overlap_outputs_Apr19_2023.csv")) %>%
  rename(ROI_2 = roi_2,
         subjID = participantID)
## Parsed with column specification:
## cols(
##   participantID = col_character(),
##   ROI_1 = col_character(),
##   roi_2 = col_character(),
##   overlap_dices_coef = col_double()
## )
MD_physics_APC_bilateral <- overlap_data_prelim %>%
  filter((
    ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" &
      ROI_2 == "MD_ROI_3_13_bilateral"
  ) |
    (
      ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" &
        ROI_2 == "MD_ROI_3_13_bilateral"
    )
  ) %>% 
  mutate(hemisphere = "bilateral")

MD_physics_APC_byhem <- overlap_data_prelim %>%
  filter((
    ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" &
      ROI_2 == "MD_ROI_13"
  ) |
    (
      ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" &
        ROI_2 == "MD_ROI_3"
    )
  ) %>%
  mutate(hemisphere = case_when(ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" ~ "right",
                                   ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" ~ "left"))

names(MD_physics_APC_bilateral)
## [1] "subjID"             "ROI_1"              "ROI_2"             
## [4] "overlap_dices_coef" "hemisphere"
names(MD_physics_APC_byhem)
## [1] "subjID"             "ROI_1"              "ROI_2"             
## [4] "overlap_dices_coef" "hemisphere"
MD_physics_APC_bilateral <- rbind(MD_physics_APC_bilateral, MD_physics_APC_byhem)
overlap_summary <- ggplot(MD_physics_APC_bilateral, aes(hemisphere, overlap_dices_coef)) +
  geom_boxplot() +
  geom_point(fun = "mean", stat = "summary", shape = 18, colour = "red", size = 4) +
  geom_jitter(height = 0, width = .1) +
  ylab("Dice's Coefficient") +
  xlab("APC ROI (Bilateral = Confirmatory)") 
  # geom_text_repel(aes(label = subjID), size = 3) +
  # ggtitle("Dice's Coefficient between LSMG and RSMG \n to bilateral APC, left APC, and right APC")
overlap_summary

APC_SMG_overlap_results <- MD_physics_APC_bilateral %>%
  group_by(hemisphere) %>%
  summarise(mean_DC = mean(overlap_dices_coef),
            median_DC = median(overlap_dices_coef),
            range_low = range(overlap_dices_coef)[1],
            range_high = range(overlap_dices_coef)[2])

MD_physics_APC_bilateral %>%
  filter(overlap_dices_coef == 0) %>%
  group_by(hemisphere) %>%
  summarise(n = n())
MD_physics_APC_bilateral %>%
  filter(overlap_dices_coef != 0,
         hemisphere == "bilateral") %>%
  group_by(subjID) 
MD_physics_APC_byhem %>%
    group_by(hemisphere) %>%
  summarise(mean_DC = mean(overlap_dices_coef),
            median_DC = median(overlap_dices_coef),
            range_low = range(overlap_dices_coef)[1],
            range_high = range(overlap_dices_coef)[2])
MD_physics_APC_byhem %>%
  filter(overlap_dices_coef == 0) %>%
  group_by(hemisphere) %>%
  summarise(n = n())
# write_rds(APC_SMG_overlap_results, path = here("outputs/univariate_results/APC_SMG_overlap_results.Rds"))

9 Part 6: Effects per task

single_betas_pertask_results <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds")) 
  # filter((domain == "both" & extracted_run_number == "all_runs") | (domain != "both" & extracted_run_number == "runs12"))
tasks <- c("infer-constraint", "agent-solidity")
single_betas_pertask_results_focal <- single_betas_pertask_results %>% 
  filter(focal_region == 1,
         effect == "event") %>%
  filter((extracted_run_number == "all_runs" & task %in% tasks) | (extracted_run_number == "runs12" & !task %in% tasks)) %>%
  select(extracted_run_number, ROI, task, domain, d, BF, B, SE, CI_95_lower, CI_95_upper, star)

single_betas_pertask_results_focal$task <- factor(single_betas_pertask_results_focal$task, levels = c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint"))

single_betas_pertask_results_focal$ROI <- factor(single_betas_pertask_results_focal$ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_bilateral", "MT_bilateral"))


pertask_b_pertask <- ggplot(single_betas_pertask_results_focal,
       aes(ROI, B, colour = task, group = task, label = star)) +
  geom_point() +
  geom_errorbar(
    aes(ymin = B - SE, ymax = B + SE, 
    colour = task),
    position = position_dodge(width = .9),
    width = .2
  ) +  
  geom_text(size = 4, aes(ROI, B + SE + .02)) +
  # geom_line() +
  geom_hline(yintercept = 0) +
  facet_wrap( ~ task, nrow = 1) +
 scale_colour_manual(values = c("#25a88b", "#19626a",  "#f24b00", "#deae1e", "#6600cd", "#fba0d0")) +
  ylab("B and SE") +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  )) +
  coord_cartesian(ylim = c(-.6, .8))

pertask_b_perROI <- ggplot(single_betas_pertask_results_focal,
       aes(task, B, colour = task, group = task, label = star)) +
  geom_point() +
  geom_errorbar(
    aes(ymin = B - SE, ymax = B + SE, 
    colour = task),
    position = position_dodge(width = .9),
    width = .2
  ) +  
  geom_text(size = 4, aes(task, B + SE + .02)) +
  # geom_line() +
  geom_hline(yintercept = 0) +
  facet_wrap( ~ ROI, nrow = 1) +
 scale_colour_manual(values = c("#25a88b", "#19626a",  "#f24b00", "#deae1e", "#6600cd", "#fba0d0")) +
  ylab("B and SE") +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  )) +
  coord_cartesian(ylim = c(-.6, .8))
pertask_b_pertask

pertask_b_perROI

10 Part 7: Region by region univar effect size

This analysis accounts for low level visual confounds.

data_for_modeling <- single_betas %>%
  filter(event != "fam",
         top_n_voxels == "100",
         domain != "psychology-environment") %>%
  filter(extracted_run_number == "run1" | extracted_run_number == "run2") %>%
  filter(!str_detect(ROI_name_final, "bilateral"))

data_for_modeling$ROI_name_final <- as.factor(data_for_modeling$ROI_name_final)

ROIs <- sort(unique(data_for_modeling$ROI_name_final))

data_for_modeling$event <- relevel(data_for_modeling$event, ref = "unexp")
modelsummaries_init <- data.frame()

for (ROI in ROIs) {
  curROI <- ROI

  # FIRST COMPUTE EVENT EFFECT SIZES, FOR EACH DOMAIN
  physics_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           domain == "physics") 
  
  psychology_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           domain == "psychology-action") 

      
  # had to remove random intercept for runs (convergence))
  physics_model <- lmer(data = physics_data,
                formula = meanbeta ~ event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
  
  psychology_model <- lmer(data = psychology_data,
                formula = meanbeta ~ event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
  
  physics_event_effect_size <- lme.dscore(physics_model, physics_data, "lme4") %>%
    select(d) %>%
    slice(1)
  
  psychology_event_effect_size <- lme.dscore(psychology_model, psychology_data, "lme4") %>%
    select(d) %>%
    slice(1)
  
  # SECOND COMPUTE DOMAIN EFFECT SIZES, FOR EACH EVENT
  expected_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           event == "exp")

  
  unexpected_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           event == "unexp")
  # had to remove random intercept for runs (convergence))
  exp_model <- lmer(data = expected_data,
                formula = meanbeta ~ domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
  
  unexp_model <- lmer(data = unexpected_data,
                formula = meanbeta ~ domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
  
  exp_domain_effect_size <-
    lme.dscore(exp_model, expected_data, "lme4") %>%
    select(d) %>%
    slice(1)
  
  unexp_domain_effect_size <-
    lme.dscore(unexp_model, unexpected_data, "lme4") %>%
    select(d) %>%
    slice(1)
  
  if (physics_data$ROI_category[1] == "MD" | physics_data$ROI_category[1] == "early_visual") {
    ROI_domain <- "general"
  } else {
    ROI_domain <- "specific"
  }
  
  modelsummary <-
    cbind(
      ROI_domain,
      physics_data$ROI_category[1],
      curROI,
      physics_event_effect_size,
      psychology_event_effect_size,
      exp_domain_effect_size,
      unexp_domain_effect_size
    ) %>%
    as.data.frame()
    
      modelsummaries_init <-
        rbind(modelsummaries_init, modelsummary)
}
colnames(modelsummaries_init) <-
  c("ROI_domain",
    "ROI_category",
    "ROI",
    "physics_event_effect_size",
    "psychology_event_effect_size",
    "exp_domain_effect_size",
    "unexp_domain_effect_size"
  )
rownames(modelsummaries_init) <- NULL

# write_rds(modelsummaries_init, here("outputs/univariate_results/univar_manyregions_viscontrols_results.Rds"))
univar_manyregions_results <- read_rds( here("outputs/univariate_results/univar_manyregions_viscontrols_results.Rds"))
plot1 <- ggplot(univar_manyregions_results, aes(psychology_event_effect_size, physics_event_effect_size)) +
    geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  coord_cartesian(xlim = c(-0.5, 1), ylim = c(-0.5, 1)) +
  geom_point(aes(colour = ROI_category), size = 3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ROI_domain, scales = "free") +
  scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
  # coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
  xlab("VOE Effect Size (d), Psychology Events") +
  ylab("VOE Effect Size (d), Physics Events") +
  ggtitle(label = "Exp 2 VOE effects for each domain across regions") +
  geom_text_repel(aes(label = ROI), size = 3) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        legend.position = "none")


plot2 <- ggplot(univar_manyregions_results, aes(exp_domain_effect_size, unexp_domain_effect_size)) +
    geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  coord_cartesian(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
  geom_point(aes(colour = ROI_category), size = 3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ROI_domain, scales = "free") +
    scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
  # coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
  xlab("Domain Effect Size (d), Expected Events") +
  ylab("Domain Effect Size (d), Unexpected Events") +
  ggtitle(label = "Exp 2 Domain effects for each event across regions") +
  geom_text_repel(aes(label = ROI), size = 3) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        legend.position = "none")
plot1
## `geom_smooth()` using formula = 'y ~ x'

plot2
## `geom_smooth()` using formula = 'y ~ x'

# instead of testing whether the linear relationship between x and y is 0, 
# we test for independence instead. H0 is that x and y are independent; F_{XY}(x,y) = F_X(x) F_Y(y).
observed_cors_ind <- univar_manyregions_results %>%
  group_by(ROI_domain) %>%
  summarise(
    cor_domain_within_event =
      np.cor.test(
        exp_domain_effect_size,
        unexp_domain_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$estimate,
    
    p_domain_within_event =
      np.cor.test(
        exp_domain_effect_size,
        unexp_domain_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$p.value,
    
    cor_event_within_domain =
      np.cor.test(
        psychology_event_effect_size,
        physics_event_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$estimate,
    
    p_event_within_domain =
      np.cor.test(
        psychology_event_effect_size,
        physics_event_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$p.value
  )
observed_cors_ind
# write_rds(observed_cors_ind, path = here("outputs/univariate_results/study2_univar_manyregions_observed_cor.Rds"))

10.0.1 Bootstrap stats

set.seed(2020)
## First define a function to work out the difference of corrs:
diff_corr <- function(data, indices) {
  data <- data[indices,]
  cor1 <-
    np.cor.test(DS_results$exp_domain_effect_size,
                DS_results$unexp_domain_effect_size,
                alternative = "two.sided",
                independent = TRUE,
                parallel = TRUE)$estimate
  cor2 <-
    np.cor.test(
      data$psychology_event_effect_size,
      data$physics_event_effect_size,
      alternative = "two.sided",
      independent = TRUE,
      parallel = TRUE)$estimate

    return(cor1 - cor2)
}

10.0.2 Domain specific

DS_results <- univar_manyregions_results %>%
  filter(ROI_domain == "specific")
doMC::registerDoMC(cores = 2)  # for running in parallel

# Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DS <- boot(
#   data = DS_results,
#   R = 4000,
#   statistic = diff_corr,
#   stype = "i"
# )

# saveRDS(res_boot_DS, here("outputs/univariate_results/study2_univariate_dsregions_4000perms_confirmatory.Rds"))

res_boot_DS <- read_rds(here("outputs/univariate_results/study2_univariate_dsregions_4000perms_confirmatory.Rds"))

## Retrieve the empirical 95% confidence interval:
ds_results <- append(boot.ci(res_boot_DS, type = "perc", conf = 0.95),boot.pval(res_boot_DS))

# saveRDS(ds_results, here("outputs/univariate_results/study2_univariate_dsregions_summary.Rds"))

10.0.3 Domain general

DG_results <- univar_manyregions_results %>%
  filter(ROI_domain == "general")
## Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DG <- boot(data = DG_results,
#                  R = 4000,
#                  statistic = diff_corr,
#                  stype = "i")

# saveRDS(res_boot_DG, here("outputs/univariate_results/study2_univariate_dgregions_4000perms_confirmatory.Rds"))
res_boot_DG <- readRDS(here("outputs/univariate_results/study2_univariate_dgregions_4000perms_confirmatory.Rds"))

plot(res_boot_DG)

## Retrieve the empirical 95% confidence interval (adjusted to 90% because of one-tailed prediction):

dg_results <- append(boot.ci(res_boot_DG, type = "perc", conf = 0.95),boot.pval(res_boot_DG))
# saveRDS(dg_results, here("outputs/univariate_results/study2_univariate_dgregions_summary.Rds"))